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Abstract 

Carbon nanotubes naturally tend to form crystals in 
the form of hexagonally packed bundles or ropes that 
should exhibit a transversely isotropic constitutive 
behavior. Although the intratube axial stiffness is on the 
order of 1 TPa due to a strong network of delocalized 
bonds, the intertube cohesive strength is orders of 
magnitude less controlled by weak, nonbonding van der 
Waals interactions. An accurate determination of the 
effective mechanical properties of nanotube bundles is 
important to assess potential structural applications 
such as reinforcement in future composite material 
systems. A direct method for calculating effective 
material constants is developed in the present study. 
The Lennard-Jones potential is used to model the 
nonbonding cohesive forces. A complete set of 
transverse moduli are obtained and compared with 
existing data. 

Introduction 

Future nanostructured composite materials are 
expected to incorporate carbon nanotube reinforcement 
either dispersed individually or as nanofilamentary 
bundles or ropes yielding unprecedented mechanical 
properties. A carbon nanotube is a cylindrical molecule 
composed of single or multiple walls of graphene 
sheets. These sheets are, in turn, composed of 
hexagonal units or graphene rings of carbon atoms that 
are bonded through highly stable sp 2 hybridized 
orbitals. A typical carbon nanotube is schematically 


depicted in Figure 1 while Figure 2 contains a cross- 
section of a bundle ensemble of individual nanotubes . 

Numerous studies have been made to analytically 
and experimentally determine the elastic properties of 
individual nanotubes. Axial Young's moduli on the 
order of 1 TPa have been measured using atomic force 
microscopy (AFM) 2 and thermal vibrations . Analytical 
studies have utilized ab initio calculations 4 , tight- 
binding methods 5 , molecular dynamic simulations 
(MD) 6 and lattice dynamics 78 . These elastic properties 
are entirely based on the strong intratube valence forces 
of the carbon-carbon bonds. 

Nanotube ensembles, however, typically form 
hexagonally packed crystal configurations in which the 



Figure 1. Single- walled nanotube. 
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intertube force interactions are due exclusively to non- 
bonding van der Waals effects which are much weaker 
than the valence forces and are highly nonlinear. Less 
consideration has been given to the transverse 
mechanical properties of nanotube bundles which 
depend on a good description of these non-bonding 
interactions. These intertube cohesive properties are of 
special interest for use in predicting the properties of 
carbon nanotube polymer composites 9 and fibers of 
woven nanotubes 10 . Selected moduli of nanotube ropes 
have been calculated with a continuum model based on 
the integrated average of the discrete Lennard-Jones 
potential 11 , MD simulation using the Tersoff-Brenner 
potential 12 and lattice dynamic methods 13,14 . 

For the present study, a direct summation of atom- 
pair potentials is developed to avoid any simplifications 
made to the nonlinear van der Waals interactions and to 
permit optimum flexibility in representing 
discontinuous or curved nanotubes. Additionally, a 
direct method permits irregular and disrupted nanotube 
lattices to be modeled to simulate damage. Because the 
fundamental constituents of nanotube bundles are only 
resolvable at nanometer length scales, analyses to 
predict macroscopic properties must necessarily merge 
concepts and techniques from continuum elasticity 
theory and discrete molecular simulation. The 
methodology developed herein combines a unit cell 
continuum model with molecular static calculations to 
determine effective moduli in aligned carbon nanotube 
bundles. 

The Lennard-Jones potential is utilized to simulate 
the van der Waals interaction forces among carbon 
atom-pairs in aligned carbon nanotube arrays. An 
achiral “zig-zag” configuration is assumed for the 
carbon nanotubes and the tube radius is assumed small 
such that the cross-section can be considered rigid. 

The objective of this work is to formulate a 
molecular mechanics model to understand and predict 
the elastic properties of carbon nanotube bundles. This 
method utilizes a unit cell approach for determining 
selected transverse moduli of a hexagonal crystal of 
aligned nanotubes. The predicted moduli are shown to 
exhibit the transverse isotropy anticipated for a material 
possessing hexagonal symmetry. The predicted moduli 
are contrasted with available published data. 

Modified unit cell formulation 

In elasticity analyses, the method of unit cells has 
been developed to determine the effective properties of 



Figure 2. Typical nanotube bundle 1 . 


heterogeneous materials by identifying and analyzing 
convenient domains of repeating microstructure. 
Applied to the determination of effective continuum 
elastic moduli of bundles of aligned carbon nanotubes, 
a minimal repeating unit of nanotubes is defined and 
subjected to continuous field deformation modes. An 
assemblage of a primary unit cell of nanotubes with 
surrounding image regions that are required in applying 
periodic boundary conditions (pbc) is presented in 
Figure 3. In this approach the nanotubes are regarded as 
sources of potential energy that require interactions 
with nanotubes outside the unit cell. Under periodic 
boundary conditions, cells of nanotubes in the 
transverse plane and nanotube segments in the axial 
dimension are treated as images of the constituents 
within the cell and used in the calculation of potential 
energy. This permits interactions between atom-pairs 
across the boundary to avoid introducing discontinuities 
in the force field. In general, these conditions ensure 
conservation of mass and energy, avoid surface or 
boundary effects, and mathematically give the primary 
unit cell a strict periodicity such that it can be 
considered to represent an infinite ensemble of 
molecules. 

The basic approach of subjecting a molecular 
ensemble to fundamental strain modes and recovering 
effective moduli from energy measures has been used 
in molecular dynamic simulations 15 ' 16 . By combining 
concepts from continuum elasticity and molecular 
dynamics, these representative units will be referred to 
herein as ‘pbc-unit cells’. Pbc-unit cells can be 
constructed of arbitrary order but, with a proper 
definition of the repeat geometry, the unit energy of the 
primitive cell remains the same. Therefore the lowest 
order cell is used for computations. A minimum-order 
hexagonal pbc-unit cell is shown in Figure 4. 
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The initial equilibrium configuration of the 
hexagonal unit cell is determined by minimizing the 
energy of the system as the nanotubes are moved 
radially outward from a fixed center. This establishes 
the equilibrium radius, R e(/ , and the nanotube center-to- 
center separation distance, S , as shown in Figure 5. 

A rigorous definition of the pbc-unit cell 
dimensions is required to ensure invariance of the unit 
energy with cell size. The required planar area of the 
pbc-unit cell is given by 


axial dimension. The total area of the enclosed 
graphene units is then equated to the surface area of a 
perfect cylinder having the same radius as the nanotube. 
This is expressed as 

NAjK^ +l)=2nR m d, lf (2) 

where N is the number of graphene units around the 
nanotube circumference, A gu is the area of an individual 
graphene unit, and R tu is the nanotube radius. The 
radius of an (N,0) nanotube can be derived as 


J3 [~7lf 1 

R =— — bSec — 1 

4 2 N 


where M is the number of nanotubes within the unit 
cell. The effective depth of the pbc-unit cell is obtained 
by first selecting a number of repeat units 
(cicumferential rings of graphene), K seg , and adding one 
additional unit to account for image segments in the 


where b is the carbon-carbon bond length. With the 
area of a graphene unit given by A gu =3>/3b~/2, 
Equation (2) yields the effective depth of the pbc-unit 
cell as 



a) (2,3)-plane section b) (l,2)-plane section 

Figure 3. Pbc-unit cell showing outside periodic image nanotubes within the (2,3)-plane 
and image sections of nanotubes in the (l,2)-plane. 



Figure 4. Miminal hexagonal pbc-unit cell Figure 5. Equilibrium radius definition 

with surrounding image cells. ^ or a hexagonal cell. 


3 

American Institute of Aeronautics and Astronautics 



d ef , = — NbCosI 

K 




1- 


j_y 

N j 


(k, + 1) 


Material constitutive relationship 


modes to the crystal, and computing the potential 
(4) energy due to atom-pair interactions as a function ot the 
deformation kinematics. Effective elastic constants are 
then determined from the variation in the system 
potential energy as 


The force field within a nanotube crystal consists 
of a combination of strong linear bonding forces acting 
within the nanotube and weak non-bonding forces 
acting between adjacent nanotubes. This disparity 
between the magnitude of interatomic forces leads to a 
highly anisotropic constitutive relation. The hexagonal 
symmetry of the nanotube crystal dictates a transversely 
isotropic stress-strain relation given by 

^11 = + C |2 E 2 2 ^12 ^33 

a 22 “ C, 2 e n + C 22 e 22 + C 23 e 33 
a 33 =C 12 £,, +C 23 e 22 +C 22 e„ (5) 

o,, — 2 C w 6 l2 

=20^6,3 

ct 23 = (C :: “ C,,)e,, 

Additional relationships between the material constants 
are given in Reference 17 for a hexagonal system as 


where Cy is the material stiffness, V is the pbc-unit cell 
volume, O is the potential energy, and £ k is an applied 
strain mode. 

Strain modes are applied to the nanotubes in the 
crystal by the imposition of specific deformation fields. 
The G 2 3 shear modulus for a hexagonally packed 
nanotube array is calculated using a pbc-unit cell 
subjected to a pure shear strain mode as shown in 
Figure 6. The magnitude of the shear strain is given by 
twice the shear angle or y 23 = 2 0 . During a progressive 
deformation with increasing 0, the potential energy is 
computed by summing all atom-pair interactions 
between adjacent nanotubes. The G23 shear modulus is 
then obtained from the elastic strain energy, O , using a 
finite difference approximation as 


K 23 = (C„ +CJ/2 
G,, = (C 22 -Cj/2 


(6) 


1 d 2 <t> _ A : <t» 

vd aTf” A^d c(1 (a, Y; ,)' 


(B) 


This material is uniquely defined by five independent 
quantities where the ‘T direction is directed along the 
tube axis. The current effort will focus on the prediction 
of normal and shear properties in the transverse plane. 

Analysis methodology 

The methodology used to determine selected 
nanotube crystal properties involves defining an 
appropriate pbc-unit cell, applying selected strain 


where A^n is the planar area of the pbc-unit cell, d eff is 
the effective length of the nanotubes, A| and A 2 are first 
and second order central difference operators, 
respectively . 

The bulk modulus is computed by applying a 
dilatational strain as shown in Figure 7. Because the 
strain in the axial dimension, e n , is assumed to be zero, 
the dilation is defined as e = £,, +S 33 with 
£^=e 33 =e. The bulk modulus is obtained by 


Q-Q.Q.Q 
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Figure 6. Imposed shear deformation on hexagonally packed nanotube array. 
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applying Equation ( 7 ) using a strain given by 2 e. The 
calculation of the E n Young’s modulus and the 
Poisson’s ratio i) 23 is performed by applying e 22 strain 
increments in the 2-direction and repositioning of the 
tubes in the 3-direction to minimize the energy. The 
transverse repositioning of the tubes perpendicular to 
the load axis directly gives a measure of the Poisson 
contraction/expansion yielding £ 2 3- This is depicted in 
Figure 8. 



Figure 7 . Imposed dilatational strain on hexagonally packed nanotube array. 



Figure 8. Calculation of Poisson’s ratio. 


der Waals radius, and r Vi is the separation distance 
between the i th and j th atoms in a pair. The term 
represents the attractive contribution to the van der 
Waals forces between neutral molecules. It includes 
permanent dipole-dipole interactions, the induction 
effect of permanent dipoles, and instantaneous dipole 
induced dipole interactions which are sometimes 
referred to as the London dispersion forces. The other 
component of the van der Waals interactions mimics 
the repulsion between overlapping electron clouds and 
is modeled by the r y 12 term which is short ranged 18 19 . 


The E 22 modulus is obtained using Equation ( 7 ) 
and the Poisson ratio is determined as u 23 = “£23/^22 • 
The C 22 stiffness coefficient is obtained by applying the 
same deformation mode defined in Figure 8, but no 
energy minimization is performed to induce a 
transverse Poisson deformation. 

Potential energy calculations 


The intertube forces are typically modeled by the 
Lennard-Jones potential to represent van der Waals 
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Results and discussion 

In the present study, a crystal of single-walled 
nanotubes is analyzed. Each nanotube has a diameter of 
0.94 nm and an achiral zig-zag conformation with 12 
graphene units around the circumference. The carbon- 

o 

carbon bond length is prescribed as 1.42 A . Nanotubes 
of this size may be considered rigid in the transverse 
direction 1213 . Therefore, the only degree of freedom 
included in the deformation kinematics is the relative 
motion of the nanotube center, and the only 
contribution to the potential energy changes with 
imposed motion is computed using the Lennard-Jones 
potential. The parameters used in the Lennard-Jones 
potential are q = 34.0 (K) and a = 0.3406 nm ' . A 
comparison between predicted elastic moduli using the 
current direct method and results obtained using 
alternate approaches is presented in Table 1. From the 
limited published results it is clear that there is a wide 
variation in predicted elastic moduli for nanotube 
bundles. 


0,, . Substituting the calculated values for the Young's 
modulus E 22 and Poisson's ratio l) 23 into the resulting 
expression, the transverse normal stiffness is obtained 
as 

F 60 3 

C = 22 -_= ou — = 68.2GPa (11) 

22 l-o; 3 1 - 0.34 2 

which closely agrees with the independently calculated 
value of Cn Next, computing C23 as u 23 C 22 = 23.2 
GPa and applying the relations given by Equations (6), 
it is found that the transverse plane stiffnesses are 
recoverable from the computed shear and bulk moduli 
as 

C22 = K,3 + G 23 = 45.8 + 22.5 = 68.3 GPa 

(13) 

C 23 = K 23 - G 23 = 45.8 - 22.5 = 23.3 GPa 

Thus, the computed elastic constants are completely 
self-consistent for a transversely isotropic material as 
required. The values reported in Reference 13 (Table 1) 


Table 1. Comparison of predicted elastic constants. 


Elastic Constant 

Direct Method 

Reference 13’ 

Reference 14** 

Reference 12 

Bulk Modulus K.2 3 (GPa) 

45.8 

42.0 

18.0 

33.6 

Shear Modulus G 2 3 (GPa) 

22.5 

5.3' 

- 

— 

Young’s Modulus E 22 (GPa) 

60.3 

17.0 

- 

-- 

Normal Stiffness C 22 (GPa) * 

68.3 

42.0 

78.0 

— 

Poisson Ratio 

0.34 

0.75 

-- 

” 


* Value derived using relationship in Equation (6). 

Results generated using (7,7) helical nanotubes with diameter — 0.94 nm. 


All the results listed in Table 1 for the present 
analysis were computed independently, none were 
derived from a subset of other values. For a 3-D solid 
exhibiting hexagonal symmetry, the expression for the 
transverse Young's modulus, E 22 » is given by 

p _ (€33 ~ C 32 )[(C 33 + C 32 )C n - 2C 31 ] 

“ (C„C 11 -C^) 

Because the axial modulus of the nanotubes in the 
crystal is generally two orders of magnitude greater 
than the transverse moduli, we may take 
C„ »C„,C,,,C„. Using the relationship 

C„ = 0,3 C„ and equating Poisson ratios under 
transverse plane isotropy as U ;3 = D J2 , Equation (10) 
reduces to a simple relationship between C22. E22 and 


are also self-consistent for transverse isotropy and 
compare favorably with the present analysis for the 
prediction of the bulk modulus, K23. In contrast, an 
analytical result presented in Reference 12 (Table 1) is 
slightly lower. The calculation of the normal stiffness, 
C 22 , which is obtained from applying a similar 
deformation as that used to compute the bulk modulus 
(applying only £22 instead of both E22 and £53 with £22 = 
£33), is distinctly different from Reference 13. The 
present analysis yields a stiffness closer to that of 
Reference 14 (Table 1) in which nanotubes of the same 
radius were used but with different helicity. The effect 
of helicity is discussed in Reference 13 in which it is 
shown for several cases that the configuration of 
graphene units on the nanotube surface has a negligible 
effect for small radius tubes ( R < 16 A ). Furthermore, 
from results presented in Reference 13, C23 rnay be 
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calculated as C 23 = \) 23 C 22 =31.5 GPa. If one applies 
the first of the relations given in Equations (6), one 
obtains a value for G 23 of 5.3 GPa. This shear modulus 
is comparable to the shear modulus associated with 
parallel planes in graphite which is experimentally 
measured as 4.0 GPa 21,22 and indicates a significant 
discrepancy between predictions using the present 
method and the lattice dynamics approach used in 
Reference 13. In the direct summation and lattice 
dynamic methods, the physics of cohesion are 
identically represented by the same parameterization of 
the Lennard-Jones potential. Possible differences in 
predictions using lattice dynamics may be due to the 
inherent integral averaging of force constants used in 
the lattice dynamical matrix and the a priori selection of 
interacting nearest-neighbor atoms used in defining the 
primitive lattice cell, both of which are avoided in the 
direct method. 

A potential source of inaccuracy effecting all methods 
is the form of the Lennard-Jones potential function 
itself. The Lennard-Jones potential was originally 
developed for noble gases and is known to produce 
poor results in other applications including graphite. It 
gives good results for the C 33 modulus (interplanar 
separation), but the C 44 parallel plane shear modulus is 
underpredicted by an order of magnitude^. Alternative 
potentials have been proposed 23 that account for the 
delocalized electronic configuration of graphite and 
yield accurate predictions for both the transverse 
normal and shear moduli. Additional study is warranted 
to assess'the spatial interactions of delocalized bonds in 
carbon nanotubes that may be underestimated using a 
spherical Lennard-Jones model. 

Conclusions 

A consistent method has been formulated for and 
applied to computing effective transverse mechanical 
properties of nanotube crystals. The method is based on 
specifying a unit cell configuration with periodic 
boundary conditions, applying a deformation field 
associated with a particular strain mode, and utilizing a 
direct summation procedure to compute changes in 
potential energy from which an effective elastic 
modulus may be obtained. The disparity between 
reported predictions of mechanical properties that 
depend exclusively on van der Waals cohesion and the 
paucity of available experimental data suggest that 
much additional investigation is warranted in this area. 
The development of a more realistic representation of 
van der Waals interactions between nanotube surfaces 
may be required to correlate analytical predictions with 
future experimental measurements of nanotube crystal 
properties. 
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